大規模データマイニングでのモデル探索手法:K-sample plot

巨大地震が日本を襲い、皆不安を感じながら生活していると思います。

そんな中せめて自分に出来ることをしようと思い、研究してきた内容をブログに記します。




サンプル数が大規模なデータニューラルネットワークとかサポートベクターマシンとかをしたくても、

時間がかかってしょうがない!ってときに参考にしてみて下さい。


近年、特にweb関係の業界ではデータデータをいくらでも記録できるようになったため、データの規模が非常に大きくなっています。

大規模データに機械学習、マシーンラーニング適用したいという要望は高まっています。

そういうときはデータからサンプリングして性能を確かめることが多いと思います。


ですがそんな時は、

  • サンプル数はどれくらいがいいの?
  • 一回やっただけじゃ真の性能は分からないよね?
  • しかもクロスバリデーションしなきゃいけないし。。

などのような事を疑問に思うでしょう。





今回紹介するK-sample plot(以下、K's plot)のアイディアはとても単純で、「サンプル数を少しずつ増やしながら性能をプロットする」というものです。




これから事例を紹介します。


あるデータがあって、結果変数yと説明変数x1とx2が100万サンプル(人)であったとします。

yとx1は線形な関係、yとx2は2次的な関係があったとし、次のように乱数発生させます。

set.seed(1)
x1   <- rnorm(1000000)
set.seed(2)
x2   <- rnorm(1000000)
set.seed(3)
y    <- 2*x1 + x2**2 + rnorm(1000000)
Data <- data.frame(x1 = x1, x2 = x2, y = y)


このデータに対して以下の3種類のモデルを仮定して適用してみます。

  1. 線形モデル  :y=x1+x2
  2. 2次線形モデル:y=x1+x2+x2^2
  3. SVMモデル  :y=svm(x1, x2)

モデル1とモデル2は線形モデルなので、100万サンプル全体を使ってもそんなに時間はかかりません。

だいたい6~7秒くらいでした。

AllModel1 <- lm(y~., data = Data)		#6.24second
summary(AllModel1)


しかし、モデル3はSVMなのでとても時間がかかります。

AllModel3 <- svm(y~., data = Data)

今30分くらい経ってますがまだ終わってませんorz



大規模データに対して不用意に機械学習を実行すると悲惨です。。。





そこで、K's plotを行ってみます。

天下り式ですが、結果からお見せします。

100万サンプルから10~1,000サンプルを抜き出し、モデル1~3を当てはめてR2乗を計算しています。

そしてR2乗とサンプル数の関係をプロットしています。

library(e1071)	#SVM
X1      <- data.frame(x1 = x1, x2 = x2)
X2      <- data.frame(x1 = x1, x2 = x2, x3 = x2**2)
y       <- y
Ksample1 <- c(seq(10,100,10), seq(100,1000,100))

KsResult1 <- KsamplePlot(X1, y, Ksample1)			#2.39second
KsResult2 <- KsamplePlot(X2, y, Ksample1)			#2.39second
KsResult3 <- KsamplePlot(X1, y, Ksample1, Method = "svm")	#2.45second

3つのモデルのK's plotの結果です。

横軸がサンプル数で、縦軸が予測性能です(plateauが推定値)。

プロットの数だけ予測を行っていますので、SVMも何十回とやっているのです。


1. 線形モデル  :y=x1+x2

f:id:isseing333:20110313120423j:image


2. 2次線形モデル:y=x1+x2+x2^2

f:id:isseing333:20110313120424j:image


3. SVMモデル  :y=svm(x1, x2)

f:id:isseing333:20110313120425j:image


どのグラフも計算に2~3秒しかかかってません!



サンプル数は10~1000の間で変化させています。10~100までは10ずつ、100~1000までは100ずつ変化させています。

予測性能は1-R二乗なので、低いほど良いモデルになります。




このプロットは検出力とサンプル数の関係を擬似的に表したものであり、サンプルが増えるとどのように予測性能が上がるかを視覚的にチェックできます

この3枚のプロットを見ると、モデル2が最も良いことが分かります。



モデル2はサンプルが200あたりで予測性能が一定になっているように思います。

しかしモデル3はまだ一定ではなく、まだ下がっている途中のようにも見えます。

そこでK-sampleの数を2000まで増やしてモデル3のK's plotを描いてみます。

Ksample2 <- c(seq(10,100,10), seq(100,1000,100), seq(1000,2000,100))
KsResult3 <- KsamplePlot(X1, y, Ksample2, Method = "svm")	#7.54second

f:id:isseing333:20110313125842j:image

予想通り予測性能が少し上がりました。

しかし、モデル2には届かないようです。



このデータの場合はSVMより2次項を入れたモデルの方が良いということが示唆できます。

もう一つ分かる事は、このデータではSVMの方が予測性能の収束が遅いという事。

予測性能がいつ収束するかというポイントは、モデル探索する時に非常に興味がある点です。

場合によっては"elbow point"の推定にも使えるのではないかと思います。





さて、ここからK's plotのプログラムを紹介します。

KsamplePlot <- function(X, y, Ksample = c(seq(40,100,10), 150, seq(200,1000,100)), 
			Method = "lm", Caret = "No", size = 5, Type = "numeric"){
  library(e1071)
  library(caret)
  library(nnet)
  library(MASS)
  library(caTools)
  library(mda)
  library(glmnet)
  library(randomForest)
  
  Rsq <- rep(NA, length(Ksample)*3)
  for(j in 1:3){
    for(i in 1:length(Ksample)){
      ksample  <- Ksample[i]
      Sample   <- sample(1:nrow(X), ksample)
      KSamplex <- X[Sample, ]
      KSampley <- y[Sample]
      
      
      #---Half CV
      Train     <- sample(1:ksample, ksample / 2)
      Trainx    <- KSamplex[Train, ]
      Trainy    <- KSampley[Train]
      TrainData <- data.frame(Trainx, y = Trainy)
      
      Varidx    <- KSamplex[-Train, ]
      Varidy    <- KSampley[-Train]
      VaridData <- data.frame(Varidx, y = Varidy)
      
      
      #------numeric
      if(Type == "numeric"){
        #---develop model
        if(Caret == "No"){
          if(Method == "lm")     Model <- lm(y~., data = TrainData); Varidx <- data.frame(Varidx)
          if(Method == "svm")    Model <- svm(y~., data = TrainData)
          if(Method == "nn")     Model <- nnet(y~., data = TrainData, size = size, linout = T)
          if(Method == "rf")     Model <- randomForest(y~., data = TrainData)
          if(Method == "mars")   Model <- mars(Trainx, Trainy)
          #if(Method == "bruto") Model <- bruto(Trainx, Trainy); Varidx <- do.call(cbind, Varidx)
          if(Method == "cart")   Model <- rpart(y~., data = TrainData)
          if(Method == "lasso"){
            cvob1    <- cv.glmnet(do.call(cbind, Trainx), Trainy)
            CvLambda <- cvob1$lambda[cvob1$cvm == min(cvob1$cvm)]
            Model    <- glmnet(do.call(cbind, Trainx), Trainy, lambda = CvLambda)
            Varidx   <- do.call(cbind, Varidx)
          }
        }
        if(Caret == "Yes"){
          modeltxt <- parse(text=paste("train(y~., data=TrainData, method= Method)", sep=""))
          Model    <- eval(modeltxt)
        }
        Err    <- as.numeric(Varidy - predict(Model, Varidx))
        MSE    <- Err%*%Err
        Var    <- (Varidy-mean(Varidy))%*%(Varidy-mean(Varidy))
        Rsq[i + length(Ksample)*(j-1)] <- (Var - MSE)/Var
      }
      #------class
      if(Type == "binary"){
        #---develop model
        if(Caret == "No"){
          if(Method == "lm"){
            Model <- lda(y~., data = TrainData)
            Post  <- predict(Model, VaridData)$posterior[, 2]
            ROC   <- caTools::colAUC(Post, Varidy)
            Rsq[i + length(Ksample)*(j-1)] <- as.numeric(ROC)
          }
          if(Method == "svm"){
            Model <- svm(y ~ ., data = TrainData, kernel="polynomial", degree=3, probability = T)
            Post  <- predict(Model, VaridData, probability = T)
            ROC   <- caTools::colAUC(Post, Varidy)
            Rsq[i + length(Ksample)*(j-1)] <- as.numeric(ROC)
          }
        }
        if(Caret == "Yes"){
        }
      }
    print(paste("Ksample", j, " ", ksample, sep = ""))
    }
  }
  Ks     <- cbind(Ksample, 1 - Rsq)
  Ks     <- Ks[(0 <= Ks[, 2] & Ks[, 2] <= 1), ]
  EstPar <- KsPlot(Ks, Method = Method, Type = Type)
  list(OneExpPar = EstPar, Ksample = cbind(Ksample, 1 - Rsq))
}






KsPlot <- function(UseData, Type = "numeric", Method){
  Par    <- data.frame(k = rep(seq(0.01, 1, 0.01), 100), plateau = rep(seq(0.01, 1, 0.01), rep(100, 100)))
  Object <- rep(NA, 10000)
  for(i in 1:10000){
    k         <- Par[i, 1]
    plateau   <- Par[i, 2]
    x         <- UseData[, 1]
    y         <- UseData[, 2]
    if(Type == "numeric") yhat <- (1   - plateau)*exp(-k*x) + plateau
    if(Type == "binary")  yhat <- (0.5 - plateau)*exp(-k*x) + plateau
    Object[i] <- sum((y - yhat)**2)
  }
  EstPar <- Par[Object == min(Object), ][1, ]
  
  if(Type == "numeric") plot(UseData, xlim = c(0, max(UseData[, 1])), lwd = 2, ylim=c(0, 1), xlab="K-Sample", ylab="1 - Explained Var", main = paste("Method = ", Method, sep = ""))
  if(Type == "binary")  plot(UseData, xlim = c(0, max(UseData[, 1])), lwd = 2, ylim=c(0, 1), xlab="K-Sample", ylab="1 - AUC", main = paste("Method = ", Method, sep = ""))
  k       <- as.numeric(EstPar[1])
  plateau <- as.numeric(EstPar[2])
  XLim    <- max(UseData[, 1])
  if(Type == "numeric") Ks.Est <- (1   - plateau)*exp(-k*c(1:XLim)) + plateau
  if(Type == "binary")  Ks.Est <- (0.5 - plateau)*exp(-k*c(1:XLim)) + plateau
  lines(c(1:XLim), Ks.Est, lwd = 2)
  text(XLim*0.7, 0.8, paste("k (slope):",       round(EstPar[1], digits=2)), adj=0)
  text(XLim*0.7, 0.7, paste("plateau  :", round(EstPar[2], digits=2)), adj=0)
  EstPar
}
  • KsPlot

サンプルデータの予測データを入れるとK's plotを描く

下のKsamplePlot関数に埋め込まれているので直接は使わない

K-sampleと予測性能に対して指数モデルを当てはめて、予測性能の上がり方をチェックする


  • KsamplePlot

K's plotの本体関数

データから、サンプリング、クロスバリデーション、グラフ描画を行う。

X:説明変数行列(2変数以上)

y:結果変数ベクトル

Ksample:サンプリング数の指定

Method:回帰方法の指定、デフォルトは線形回帰、今はlmとsvmのみ実装



Medhodの拡張は簡単で、「develop model」のところの「if(Method == "lm") Model <- lm(y~., data = TrainData)」を追加するだけです。

下の行のpredictに繋げれればいいので、ランダムフォレストも追加できます。

クロスバリデーションはHalf-CV(2-fold CVの片方だけ)を3回行っています。


取り急ぎ完成しているのはここまでです。

まだデータによってはたまにエラーが出たりするので完璧じゃないのですが、使ってみて下さい。




この手順で予測モデルの構築や選択を行えば、大規模データでも素早く探索できると思います。

現在プログラムにしているのは連続値の予測だけですが、2値予測の場合はAUCを指標にすることで同じロジックを使えます


説明変数行列と結果変数ベクトルさえあればKsamplePlot関数が使えますので、是非使ってみて下さい。


これ、実は今年度出したD論の一部なんです・・……(-。-) ボソッ


また拡張したら記事にします。

それでは。



修正1

  • KsamplePlot関数の中にミスがあったので修正。
  • 連続値のオプションニューラルネットを追加
  • パッケージcaretを使えるように修正
    • caret = "Yes"、Model = "使用したモデル名"で利用できる
  • 2値予測も可能
    • Type = "binary"と指定する
    • デフォルトは線形判別分析
    • 現在SVMのみ選択可能
    • これもcaretはどうしたらいい?
    • Epi::ROCでAUCを計算しているのでグラフが一杯出てしまう

修正2

  • Epi::ROCをcaTools::colAUCに修正
    • sampleのサンプリングによってエラーが出る
    • エラーは2種類
      • yが2グループ共少なくとも1つはサンプルされてないといけない
      • svd(X, nu = 0)のエラーは不明

修正3

  • MARSの追加(brutoは動かすとRがクラッシュするので追加できず)
  • LASSOの追加(CVして最適なラムダを推定する、Ksampleは最低でも40は必要)
  • ランダムフォレストの追加

動作報告があったら教えて頂けると幸いです。


サンプルプログラム

set.seed(1)
x1   <- rnorm(1000000)
set.seed(2)
x2   <- rnorm(1000000)
set.seed(3)
y    <- 2*x1 + x2**2 + rnorm(1000000)


X1      <- data.frame(x1 = x1, x2 = x2)
X2      <- data.frame(x1 = x1, x2 = x2, x3 = x2**2)
y       <- y

KsResult1 <- KsamplePlot(X1, y)
KsResult2 <- KsamplePlot(X2, y)

ページTOPへ